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We report on the first calculations of fully relativistic bi- 
nary circular orbits to span a range of separation distances 
from the innermost stable circular orbit (ISCO), deeply inside 
the strong field regime, to a distance (~ 200 km) where the 
system is accurately described by Newtonian dynamics. We 
consider a binary system composed of two identical corotating 
neutron stars, with 1.43 Mq gravitational mass each in isola- 
tion. Using a conformally flat spatial metric we find solutions 
to the initial value equations that correspond to semi-stable 
circular orbits. At large distance, our numerical results agree 
exceedingly well with the Newtonian limit. We also present a 
self consistent determination of the ISCO for different stellar 
masses. 

PACS Numbers: 04.25.Dm, 04.40.Dg, 95.30.sf, 97.80.Fk 



Neutron-star binary systems are currently of interest as 
sources of gravitational radiation HII. The development 
of a new generation of laser interferometric and cryogenic 
gravitational wave detectors ^ has renewed interest in 
theoretical models for generating gravitational radiation. 
Binary systems composed by neutron stars and/or black 
holes are the most promising sources of detectable gravi- 
tational radiation. However, the expected signal-to-noise 
ratio is so low that the extraction of information will be 
difficult without good theoretical waveforms. Thus, it is 
important to accurately model such systems. 

At the frequencies of interest for laser interferometric 
(10-100 Hz) and cryogenic (~ 1000 Hz) detectors the or- 
bits enter the strong gravity regime. Expansion methods 
applicable in the Newtonian limit begin to breakdown 
before the strong field regime. Hence, it becomes im- 
perative to model the binary evolution with a method 
which accurately extends from the Newtonian limit to 
the strong gravity regime. In this paper we present the 
first such calculation. These calculations are for one par- 
ticular scenario, namely that of two equal-mass stars con- 
strained to rigid corotation in circular orbits. Neverthe- 
less, this represents a plausible benchmark for the ex- 
pected gravity wave signal. 

Many groups have performed numerical simulations us- 
ing different approximations to this problem. Prelimi- 
nary results have been achieved using Newtonian dynam- 
ics or post-Newtonian expansion techniques || . While 
these approximations work very well when the stars are 
very far from each other, they become unreliable when 
the distance between stars reduces to a few stellar radii. 

Recently, a new fully relativistic approximation has 



been used |Ej-|8|] to model binary neutron-star systems. 
This approximation is mainly based upon a restriction 
that the spatial part of the metric tensor is forced to be 
conformally flat (CFC). Work by Wilson and Mathews 
H and Wilson, Mathews and Marronetti reported the 
first fully relativistic calculations of quasi-stable circular 
orbits employing this approximation. The most contro- 
versial U of their results predicts that, for certain con- 
ditions, the stars could collapse into black holes prior to 
the merger. This controversy is irrelevant to the present 
work in that it has been demonstrated [8|JlO[| that this 
effect does not occur when rigid corotation is imposed. 

Baumgarte et al. Jl(J recently developed a method for 
using the CFC to compute rigidly corotating stars. They 
found circular orbits for very close separation distances 
(a few stellar radii) and estimated the point of secular 
instability of the stars. In the present work however, we 
implement a scheme to directly determine the location 
of the ISCO for systems with different masses. We also 
present numerical solutions for quasi-stable circular or- 
bits for corotating neutron stars that cover a wide range 
of separation distances between the stars. Hence for the 
first time we present the much needed calculation con- 
necting the strong field regime at the innermost stable 
circular orbit (ISCO) with the weak field region where 
the system is well described by Newtonian dynamics. 

A full discussion of the CFC method can be found in 
0. We use the (3+1) spacetime slicing as defined in 
the ADM formalism pd] , ^2[ . Utilizing Cartesian x,y,z 
isotropic coordinates, proper distance is expressed 

ds 2 = -{a 2 - I3 n (3 n )dt 2 + 2j3 n dx n dt + lns dx n dx s , (1) 

where the lapse function a describes the differential lapse 
of proper time between two hypersurfaces. The quantity 
f3i is the shift vector denoting the shift in space-like co- 
ordinates between hypersurfaces and 7^ is the spatial 
three-metric. The Latin indices go from 1 to 3. 

Using York's (3+1) formalism |l2|, the initial value 
equations can be written as follows; the Hamiltonian con- 
straint equation can be written, 



R = 16np 
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where R is the Ricci scalar curvature, K ns is the extrinsic 
curvature, and p is the mass-energy density. 
The momentum constraints have the form 11131 



D n (K m -j m K) = 87^ 



(3) 
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where D n is the three-space covariant derivative and S z 
are the spatial components of the four-momentum den- 
sity. 

The CFC method restricts the spatial metric 7^ to the 
form 



(4) 



where the conformal factor (f> is a positive scalar function. 
This approximation simplifies greatly the equations. It 
is motivated in part by the fact that the gravitational ra- 
diation in most systems studied so far is small compared 
to the total gravitational mass. The CFC is a frequently 
employed approach to the initial value problem in nu- 
merical relativity. Its application here is consistent with 
the quasi-equilibrium orbit approximation. Further jus- 
tification for its use can be found in refs. Q and fl4|| . 

The CFC leads to a set of elliptic equations for the 
metric components. Using Eq. (j^) in combination with 
the maximal slicing condition tr(K) = 0, we get the fol- 
lowing equations for <f> and (a</>) 
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V 2 (a</>) = 47T/92 
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where the V, represent flat-space derivatives and the 
source terms are 



P1 = Y 
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where po is the rest-mass density, e the internal energy 
per unit of rest mass, T the adiabatic index, and W a 
generalization of the special relativistic 7 factor Jjj. A 
solution of equation (|^) determines the lapse function 
after equation (|^) is used to determine the conformal 
factor. 

The shift vector /3 l can be decomposed ]15[: 



0i = B i --V i x ■ 



(9) 



This is introduced into Eq. (|| 
two elliptic equations 

V 2 X = V„B' 



to obtain the following 



(10) 



V z B l = 2VJn(a(j)- b )K m - 16Tra<j) 4 S i 



(11) 



An equation for the extrinsic curvature K lJ is derived 
M using the time evolution equation and the maximal 
slicing condition 



K* = %-(y i p+V i p i --8i j Vn0 n ) , (12) 
la & 

where K ij = 4> w K l i. 

We assume that the matter behaves like a perfect fluid 
with a stress-energy tensor 

r = (p„(i + £ ) + p)«v + i^ , (13) 

and use a polytropic equation of state (EOS) 

P = kp r Q , (14) 

with P the pressure and fc a constant. The results pre- 
sented here are for stars with T = 2. 

Following Baumgarte et al. |ITJ], for rigidly corotating 
stars, the fluid four velocity can be taken as proportional 
proportional to a Killing vector 
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where u> is the orbital angular frequency and $ the az- 
imuth coordinate. In this case the steady-state limit of 
the hydrodynamics momentum equation Q yields the 
relativistic Bernoulli equation [n0[ 
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q = 



l + C 



1 + n \a(l - v 2 ) 1 / 2 
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(16) 



where q = P/po, C is a constant of integration and v 
is the matter proper velocity |l6[ ]. From Eqs. ( |l3| ) and 
( p6| ) we find expressions for the proper baryonic matter 
density po and the material momentum density S l as 
functions of the fields. We choose a, (3 l , (f>, and q as the 
independent variables of our set of equations. 

The set of equations is solved numerically using an 
iterative algorithm based upon a specially designed el- 
liptic solver. This method consists of a combination 
of multigrid algorithms and domain decomposition tech- 
niques [17]] . The three dimensional spatial volume where 
the equations are solved is divided into concentric lay- 
ers. These layers are centered around the star and the 
grid resolution decreases with the distance from the stars. 
To avoid bias, each iteration cycle starts assuming zero 
angular frequency and spherical stars. The equations 
are solved iteratively until numerical convergence for the 
fields and the frequency is achieved. Our calculations 
utilize 40 grid points on average across the stellar diam- 
eter. This is about twice as large as the number used in 
previous works pQ] . This efficient use of computer mem- 
ory permits us to describe adequately the interior of the 
stars even when they are so far apart that they approach 
the regime of Newtonian point masses. 

The boundary conditions are estimated from the first 
terms in a multipole expansion of the fields. This method 
is quite accurate when the boundary surfaces are very 
far from the stars, which is the case in our computations. 
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Thus, these methods allow us to connect between the 
strong gravitational field regime close to the ISCO and 
the Keplerian regime. Domain decomposition techniques 
also provide a natural way of code parallelization, which 
reduces the processing time enormously jl7| . 

Solutions are obtained for specific values of two free pa- 
rameters, namely, the coordinate distance between stars 
and the stellar baryonic mass. The reason to choose the 
latter as a free parameter is that the baryonic mass of the 
star remains constant during the inspiral that describes 
most of the time evolution of the system. This differs 
from the approach of [Q who used central density and 
separation distances as free parameters and interpolated 
to find the solutions of constant baryonic mass. In our 
calculations we iterate to obtain a given baryonic mass. 
This avoids the numerical error associated with interpo- 
lation. Thus, we were able to construct constant baryonic 
mass sequences of orbits with a minimum number of code 
runs. 

Using the CFC approximation we found solutions to 
the initial value equations for semi-stable circular orbits 
for a binary system of identical neutron stars with 1.55 
Mq baryonic mass and 1.43 M© gravitational mass in 
isolation. The former values represent a typical neutron 
star and were obtained fixing by the value of the poly- 
tropic constant (in polytropic models physical quantities 
can be rescaled with the polytropic constant JToj]). 
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FIG. 1. Orbital frequency as a function of the coordinate 
separation between star centers. The solid line represents 
Kepler's law, the circles are the orbits corresponding to a 
sequence of stellar baryonic mass of 1.55 Mq, and the squares 
are results from Baumgarte et al. [10]. 

Figure |l| shows the angular frequency as seen by a dis- 
tant observer as a function of the coordinate distance 
between the centers of the stars (circles). The coordi- 
nate distance is not a gauge invariant quantity, but since 
it converges asymptotically to flat space separations, it 
allows us to compare our results with Kepler's law (solid 



line). We can also compare our results with those from 
Baumgarte et al. |hJ (squares) since we use the same co- 
ordinates. They estimate the onset of secular instability 
after which the stars can no longer maintain rigid coro- 
tation. They assume that the ISCO is close to this point 
since the secular instability is expected to occur just be- 
fore the dynamical instability that defines the ISCO [jl8| . 
Note that their calculations include numerical solutions 
for orbits that are inside the ISCO which are used to 
determine the turning point in the binding energy |10| . 
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FIG. 2. Binding energy of the system as a function of the 
orbital frequency as seen by a distant observer. The points 
represent the orbits of the high resolution runs and the solid 
line corresponds to a polynomial fit. We identify the lowest 
point as the ISCO of the system. 

The numerical method developed for this work shows 
the absence of solutions for circular orbits when the stars 
are too close. This situation is consistent with the exis- 
tence of a dynamical instability which defines an inner- 
most stable circular orbit for the binary system. How- 
ever, there is still the possibility that this absence is due 
to a problem in the convergence of the numerical scheme. 
The ISCO determined here is thus the stable orbit with 
the highest angular frequency (in the sequence shown in 
Fig. [y, it is at uj = 3576 rad/sec). 

The value of angular momentum at the ISCO is J — 
1.61 x 10 cm r, which gives a value of J/M G 2 = 0.94. 
Thus, the stars could merge into a Kerr black hole with- 
out the further loss of angular momentum. 

Figure ^ shows the binding energy as a function of the 
angular frequency of the system as seen by a distant ob- 
server. The binding energy is represented as one half 
the ADM mass of the binary minus the gravitational en- 
ergy of an isolated star, divided by the baryonic mass 
of the isolated star. The points correspond to the orbits 
obtained using high resolution (more than 30 grid points 
across the stellar diameter) and the solid line corresponds 
to a polynomial fit. The lowest plotted value of the bind- 
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ing energy on Fig. || coincides with the orbit we identify 
as the ISCO. In this calculation no turning point was 
observed (within the numerical error). For large separa- 
tions (low angular frequencies) the value of the binding 
energy approaches zero as expected. 

Figure H shows the ISCO angular frequency for systems 
of stars with different gravitational mass in isolation. The 
angular frequencies shown in Fig. |^ are somewhat lower 
(~ 10%) than the estimations corresponding to the sec- 
ular instability points from Baumgarte et al. JTof ] - This 
is most likely due to effects of grid resolution. 

The excellent agreement between our results at large 
separations and the values obtained from Newtonian dy- 
namics provides an additional check on the numerical 
code. 
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FIG. 3. Orbital angular frequency of the ISCO for stars 
with different gravitational masses (in units of solar masses) 
in isolation. 

We observe a slight decrease in the central density 
when the stars approach each other similar to the results 
reported in JToj ] . This is contrary to what is observed in 
fully hydrodynamical simulations. In reference Q how- 
ever it has been shown analytically and numerically that 
no increase in central density occurs for stars in rigid 
corotation. This also confirms that the effect is not an 
artifact of the CFC method. 

An analysis of the emission of gravitational radiation 
along sequences of constant baryonic mass will be the 
subject of a forthcoming article. 
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